Use of Boundary-Driven Nonequilibrium Molecular Dynamics for Determining Transport Diffusivities of Multicomponent Mixtures in Nanoporous Materials

The boundary-driven molecular modeling strategy to evaluate mass transport coefficients of fluids in nanoconfined media is revisited and expanded to multicomponent mixtures. The method requires setting up a simulation with bulk fluid reservoirs upstream and downstream of a porous media. A fluid flow is induced by applying an external force at the periodic boundary between the upstream and downstream reservoirs. The relationship between the resulting flow and the density gradient of the adsorbed fluid at the entrance/exit of the porous media provides for a direct path for the calculation of the transport diffusivities. It is shown how the transport diffusivities found this way relate to the collective, Onsager, and self-diffusion coefficients, typically used in other contexts to describe fluid transport in porous media. Examples are provided by calculating the diffusion coefficients of a Lennard-Jones (LJ) fluid and mixtures of differently sized LJ particles in slit pores, a realistic model of methane in carbon-based slit pores, and binary mixtures of methane with hypothetical counterparts having different attractions to the solid. The method is seen to be robust and particularly suited for the study of study of transport of dense fluids and liquids in nanoconfined media.


■ INTRODUCTION
Understanding and modeling of mass transport in porous media is essential in nearly all branches of natural sciences and engineering given the ubiquity of fluid flow in natural and anthropogenic solids. Particularly in engineering applications, a significant number of processes exploit the relative differences in mass transport to separate different components of a fluid mixture. By employing porous matrices, such as membranes, the optimal design of separation units can lead to significant reductions in cost and energy consumption. An early example of separation of fluids using porous media is water desalination using reverse osmosis membranes, where salt is removed from water without the need for energy-intensive distillation units. 1 Since the turn of the century, advances in technology and a perennially increasing control over nanostructure design has led to the production of porous materials with remarkable adsorption and transport properties, sometimes counterintuitive or unexpected. Fluids flowing within carbon nanotubes (CNTs) 2−5 exhibit extremely large fluxes. Structured solids such as zeolites 6,7 and metal−organic frameworks (MOFs), 8−10 with very large surface area to volume ratios, have the potential of being selective to certain molecules due to a combination of steric and energetic effects, which makes them ideal for separation and catalysis processes. Extremely thin yet strong polymer membranes have been designed to be used in nanofiltration 11,12 and more recently polymer membranes have been used as an alternative to distillation in the fractionation of crude oil. 13 The issues associated with the production of oil through unconventional tight nanoporous shale rocks is another example of the unforeseen behavior drawn by the flow of fluids through ultraconfined media. 14 Optimal design in separation processes has an enormous impact on overcoming current scientific and environmental challenges. 15 In many instances, given the nanoscale porosities of these materials, understanding the precise mechanism of mass transport is challenging. It is not known to what extent traditional kinetic and phenomenological models can be used to predict transport in these materials. Generally speaking, there are three theories used to characterize transport in porous materials. The first is based on Fick's law, which employs an empirical constant to relate the transport (mass flux) to the macroscopic density (or pressure) gradient that drives it. 16 The strategy is conceptually straightforward and commonly implemented in the analysis of experimental results.
The resulting coefficients suffer from many drawbacks, namely, a lack of transferability, and a nontrivial pressure, temperature, and concentration dependence. The Onsager formulation, based on irreversible thermodynamics, 17 recognizes that the driving force for transport is actually a gradient of chemical potential. The transport coefficients thus generated are fundamentally robust and better behaved than the Fickian counterparts, but they are conceptually challenging as chemical potentials cannot be directly measured. A final formulation, defined as the generalized Maxwell−Stefan theory, 18 suggests the description in terms of so-called corrected diffusities. There is a formal link between all these formulations. 19,20 However, in all of these descriptions, there is a need to experimentally determine the transport coefficients, i.e., there are no currently available fully predictive methods.
Experimental determination of mass transport and diffusion in nanoporous materials is challenging. 21,22 A fundamental predictive theory behind transport in nanoconfined spaces is lagging behind the bulk−fluid counterpart, mainly due to the complexities of incorporating the surface−fluid effects, which are, in many cases, dominant. 23 Thus, as a complementary tool to experiment and theoretical modeling, molecular simulations can provide useful insight into the mechanism of transport and separation in confined media. 23 Molecular dynamics (MD), Monte Carlo (MC), and a combination of both are commonly used in literature to study transport of pure fluids and multicomponent mixtures. 24−26 Calculations can be based on the analysis of systems in or away from equilibrium. 27 In the equilibrium molecular dynamics (EqMD) method, the mean squared displacement, or velocities, of individual particles and their centers of mass are used to measure the motion of the fluid. To calculate transport diffusivities, additional information, namely, adsorption isotherms, are included. Otherwise, the tracking of particle positions or velocities is employed to predict transport diffusivities of multicomponent mixtures. 18,28−32 A particularly relevant scenario to this manuscript is the calculation of transport diffusivities in nanopores. 33−36 Calculating transport diffusivities using equilibrium molecular dynamics requires very long simulation runs, 37 and additional simulations (typically Grand Canonical MC) to calculate the Darken correction factors from adsorption isotherms (∂μ/∂(ln ρ)) are required. As an alternative, nonequilibrium molecular dynamics (NEMD) techniques have been developed measuring transport by inducing a flux by external fields or by artificial gradients. 27,38 One of the first reported nonequilibrium methods is the external-field NEMD (EF-NEMD) method of Evans and Morris, 39 where an external field exerts a force on all fluid particles, generating a steadystate nonequilibrium flux. The ratio of the flux to the force is identified as the inherent transport coefficient or Onsager coefficient. This method has been successfully used to study pure systems 27,38 and binary mixtures. 40−42 A similar method of using walls to push fluid particles through pores in a pistonlike fashion has also been proposed by Wang et al. 43,44 To report transport diffusivities using these techniques, adsorption isotherms are required to calculate the Darken factors (see the "Methodology" section).
Other NEMD methods have been developed that do not require adsorption isotherms to measure transport diffusivities. These methods rely on having density gradients across the pore, with the ratio of the induced flux to the density gradient being the transport diffusivity. One of the earliest of such direct measurement methods is the gradient relaxation molecular dynamics (GRMD) 38 method where a system is initially set up with a density gradient and is then allowed to relax until fully equilibrated. The evolution of the fluid motion to an equilibrated state is then modeled by the diffusion partial different equation to measure the transport diffusivity. Given that the system never reaches steady state until fully equilibrated, it can be difficult to measure the transport diffusivity using ever changing density gradients and fluxes, leading to poor statistics. More efficient algorithms have been introduced by ensuring a steady-state gradient, and thus a steady-state flux, reducing the error in measured variables. One of the most common nonequilibrium steady-state methods is the dual control volume grand canonical molecular dynamics (DCV-GCMD) that has been used to measure transport diffusivities of pure fluids 27,45−48 and mixtures. 49,50 DCV-GCMD combines insertion and deletion methods of grand canonical MC with MD so that particles can be deleted downstream of the pore and inserted upstream to create a steady-state chemical potential and density gradient across the pore. The combination of MC with MD has clear advantages, as a steady-state gradient can be achieved and the pressures upstream and downstream can be clearly defined. However, the use of both MC and MD moves can be cumbersome with particle insertions being particularly difficult for high pressures and dense fluids. A more recent method fixing chemical potentials across a pore is that of Ozcan et al. 51 Other NEMD techniques exist that do not rely on chemical potentials, and subsequently on insertions and deletions. One relevant approach is that of Li et al., 52 where a pressure gradient is achieved using a partially reflecting membrane. A membrane is positioned in the simulation box where particles crossing the membrane in a certain direction cross it without hindrance, yet particles attempting to cross the membrane in the opposite direction have a probability of being reflected back. This leads to a pressure gradient and a steady-state flow.
Finally, the boundary-driven NEMD (BD-NEMD) method of Frentrup et al. 37 is a unique method of measuring transport diffusivities where MD is solely used to induce steady-state density gradients across nanopores. This is done by applying an external field only to a small region in the simulation box, positioned far from the pore. This tool has been extensively used in literature for pure components 37,53−57 and mixtures. 58,59 Although the BD-NEMD method has been extensively used in literature, it is not known to what extent transport diffusivities measured in this technique agree with other proven methods, such as the EqMD method. Moreover, there are significant assumptions that need to be addressed to ensure that transport diffusivities measured from the BD-NEMD and EqMD methods agree under all conditions. Moreover, although this method has been used to study multicomponent mixtures, 58 only the self-transport coefficients have been measured, and the effect of the cross-species transport coefficients have not been studied. The values of the crossspecies transport coefficients have not been previously measured using the BD-NEMD method, and it is not known to what extent transport diffusivities of multicomponent mixtures measured using this method agree with those from EqMD simulation. This work aims to address these issues.

■ METHODOLOGY
Several excellent reviews discuss the relationship between the different transport coefficients which may be directly or indirectly measured in molecular simulations. The reader is referred to them for extension on this topic. 19,60 For completeness, we will briefly discuss the most relevant expressions and relationships.
Transport under Equilibrium: Benchmark Case. In molecular simulations, a common practice is to measure the motion of molecules using the self-diffusivity of individually tagged particles of type l, which can be calculated using the Einstein relation: 61 where r l,i is the vector describing the position of the ith particle of molecular type l, N l is the total number of molecules of type l, and d is the dimension of the vector r l,i . In bulk, d has a value of 3 as particles can move in all three dimensions. However, inside porous materials, the value of d can range from 1 to 3, depending on the dimensions where particles can freely move. The angle brackets refer to an ensemble average. Although self-diffusivity is informative in describing the motion of individual particles, it does not describe the collective mobility of molecules, i.e., how the collective (center of mass) motion of molecules of type l is related to the collective motion of all molecules of type m. Given that flow in porous materials is a consequence of the collective motion of molecules, selfdiffusivity cannot be used to calculate transport of fluids apart from extremely dilute systems where particles do not experience strong intermolecular forces and thus move independently. 62 A more general equation is thus used to describe the transport of particles of type l, influenced by the collective motion of particles of type m, denoted by Λ lm : 35 For a pure component system, Λ 11 corresponds to the Maxwell−Stefan (MS) diffusivity, D̵ MS , also known as the collective, or Darken corrected, diffusivity, D c . Thus, eq 2 can be simplified for a pure component system: For multicomponent mixtures, elements of the matrix [Λ] are related, but they are not identical to the exchange coefficients of the MS diffusivity. Λ ij provides an indirect route to calculating transport diffusivities of multicomponent mixtures (eq 11). In addition to the Einstein relation, Λ ij can also be measured using Green−Kubo relations which use velocities instead of displacements. 29 Although very commonly used, calculating collective diffusivities from equilibrium molecular dynamics simulations requires many independent simulations (or one very long simulation) to get acceptable statistics. 33,62,64 This renders this method laborious and computationally inefficient.
Moreover, in order to compute transport diffusivities from equilibrium simulations, additional information is required in the form of the adsorption isotherms for the fluid in the porous material. While the equilibrium methods are accepted as the de facto "gold standard", 65 they are computationally inefficient. It is in this space that nonequilibrium models can be an alternative to calculate transport diffusivities.
In the next section, we will discuss an implementation of a boundary-driven nonequilibrium method to binary mixtures and the relationship between the transport coefficients obtained as compared to those derived from other routes.
Transport under Non-Equilibrium. Transport diffusivity, D t , is the mass transport coefficient used in the continuity and Fick's equations, relating the flux of a given species, i, in a porous medium to its concentration gradient: 36,40,66 Here, J is the n × 1 vector of fluxes of n components inside the pores, [D t ] is the n × n Fick or transport diffusion matrix, and ∇ρ is the n × 1 vector of density gradients. It is important to note that fluxes need to be measured relative to a specified frame of reference which is the porous solid and is assumed to be stationary. Moreover, it is also important to mention that ∇ρ corresponds to the density gradient of the fluid inside the porous media, i.e., the density of the adsorbate.
For a binary mixture where the flux is measured in one direction, z, the above equation can be written as 67 where D ii t is the self-transport diffusivity, i.e., the contribution to the flux of species i due to its own concentration gradient, and D ij t is the mutual diffusivity, corresponding to the contribution to the flux of species i due to the concentration gradient of species j. For a binary mixture, none of the elements of the transport diffusivity matrix are necessarily similar to each other, i.e., D ij t ≠ D ji t . Crucially, a concentration gradient is not the only source of mass flux, e.g., the Soret effect 68 describes the flux of mass driven by a thermal gradient which has been studied in bulk 68,69 and under confinement 70−72 using simulations. Similarly, consider a system in vapor−liquid equilibrium where there is a clear gradient in density; however, there is no net mass flux across the interface. As an alternative to the incongruities of the Fickian formulation, Onsager's treatment provides a fundamental starting point relating the fluxes to the underlying mass transport driving force. In an isothermal case, one could relate the fluxes to the gradients of the chemical potential difference of each species, ∇μ: 17 6) Here, [L] is the symmetric matrix of Onsager coefficients (phenomenological coefficients). An interesting inference of the Onsager treatment is that the flux of species i becomes dependent on the chemical potential gradient of both components i and j. There is an explicit recognition within this formulation that there are cross-component effects, i.e., that the flux of one component may have an impact on the The Journal of Physical Chemistry B pubs.acs.org/JPCB Article flow characteristics of the other. In particular, the cross coefficient terms become important because the density of the system becomes liquidlike and correlation effects are strong. 73 For a binary mixture which is not true for mutual transport diffusivities. The matrix Comparing eq 5 with eqs 7 and 10, the relationship between the diffusion coefficients may be expressed as where #Comps is the number of components. The term βρ i ∂μ i /∂ρ j is commonly known as the Darken factor, Γ ij , and for a pure system, it is the proportionality constant relating transport to the collective diffusivity: The Darken factor can be calculated from adsorption isotherms obtained from simulations or experimental data. A critical review of some of the alternative techniques for calculating diffusion coefficients can be found in literature. 17,27 Of particular note is the Maxwell−Stefan formalism. 18,74−76 These models are all "equivalent" and the relationship between them has been presented elsewhere. 16,77 Boundary Driven Non-Equilibrium Molecular Dynamics. Given a density gradient inside a pore and an appreciable flux, it is possible to directly assess the transport diffusivity of a fluid flowing inside a pore. To induce such gradient, we revisit the proposal of Frentrup et al. 37 as applied to pure fluids, modifying key assumptions, and extending it to multicomponent mixtures. First, a simulation box composed of two bulk reservoirs in contact with the pore is set up. Fluid particles are added until a target global density is reached and the system is left to equilibrate. An external field is then applied in a small section of the simulation box, which in this work is 2 nm wide, far enough from the pore. This region is also the boundary of the two reservoirs. This external field exerts a directional force on all particles within the small region, The Journal of Physical Chemistry B pubs.acs.org/JPCB Article pushing particles at the boundary of the two reservoirs to one side, thus increasing the concentration of the fluid in one reservoir, and depleting the amount of fluid in the other. The reason the external force is applied far away from the pore is to minimize the effects of the applied force on the flow across the pore, to ensure that the flow is being induced by the density gradient across the pore. The system eventually reaches steady state, and a steady-state concentration gradient occurs within the pore. This is illustrated in Figure 1 for a binary mixture. In terms of implementation, the directional force is the result of applying an acceleration to the particles. In this work, the same acceleration is applied to all species. Interestingly, in this method there are two bulk regions on each side of the pore. For small applied forces, the density of each bulk region remains fairly constant, and the observed density gradient is only within the pore. Applying different external forces to the boundary leads to different bulk compositions and different density gradients across pores, as can be seen in Figure 2, for a pure fluid flowing through a slit pore.
In Figure 2 (left), different forces are applied to the fluid and a density gradient develops within the pore. It can be clearly observed that the density in the adsorbed phase and in the bulk are very different. This implies that density gradients across the pore, and consequently the transport diffusivities, are different dependent on which density (bulk or adsorbed) is used. It is important to mention that commonly in literature using BD-NEMD transport coefficients are calculated using reservoir densities 37,54,56,78,79 which is only rigorous if the adsorbed densities are the same as bulk densities, i.e., for nonadsorbing systems. As will be shown later, an incorrect choice leads to significant discrepancy between measurements of transport diffusivities using the EqMD and NEMD methods.
Therefore, a key modification of the BD-NEMD method implemented in this work is to use density gradients inside the pore instead of calculating the density gradients using bulk (reservoir) concentrations. This can be particularly challenging, as the statistics inside the pore can be poor. In order to estimate the density gradients more robustly, an assumption is made that inside the pore (pore entrance and exit), the boundaries are at local equilibrium with the adjacent bulk reservoirs; thus, the amount adsorbed at the boundaries can be calculated from equilibrium simulations of a bulk region. Knowing the bulk densities on each side, the adsorbed density at the boundary inside the pore can be estimated, and the density gradient is calculated. An example of the such relationship between bulk and adsorbed density can be seen in Figure 2 (right).
Once the density gradient is calculated, the flux of each component, J i , is calculated in the middle of the pore: where A xy is the area of the plane in the pore perpendicular to the direction of the flow, t run is the total simulation running time, and N i + and N i − are the total number of particles of species i that have passed the middle of the pore in the same and opposite direction to the direction of the external force, respectively. The choice of the middle of the pore is arbitrary if the cross-sectional area of the pore does not change.
By running simulations with different forces, one could obtain different density gradients and fluxes which can be used to improve the statistics. Transport diffusivities can be assessed using the following equation: For an m component mixture, D t is the transport diffusivity matrix defined in eq 5, and ∂ρ/∂z and J are n × m matrices of concentration gradients and fluxes, respectively. n is the number of in silico experiments carried out, each with a different boundary force, and the superscript "T" denotes a matrix transpose.
One of the benefits of using the BD-NEMD technique over other nonequilibrium approaches is the fact that the force itself is not used to evaluate transport coefficients. It is applied in a region so far away from the pore that it does not affect the transport inside the pores. However, what the force does is to help in building up fluid on one side of the pore. Transport is a thus a consequence of the concentration gradient across the pore, and not the applied force.
Thus, it does not strictly matter what force is applied to either species, as long as there are measurable concentration gradients and fluxes of all species across the pore. If fluxes are only functions of concentration gradients, and not applied forces, then transport coefficients are independent of those concentration gradients and thus different forces can be applied to any species. This is in contrast with other NEMD techniques, such as the external force NEMD (EF-NEMD) method, where a force is applied to all species. In that case, the applied force is directly responsible for the transport of both species and transport coefficients are calculated by relating fluxes to applied forces. In that case, if molecules have different masses, then one might observe buoyancy effects and artifacts that might affect the measurement of transport coefficients. To overcome that particles having different masses, accelerations should be modified to ensure the same force applied to different species.
Molecular Interactions. All fluid molecules are modeled as single spheres, with the Mie potential describing intermolecular interactions between particles. Cross interactions are resolved by using the Lafitte et al. combination rules. 80 Details are provided in the Supporting Information.
Modeling Adsorption. As highlighted previously in Figure  2 (right), in the BD-NEMD method it is important to relate   61,81 and EqMD simulations have been used. The GCMC method models the pore and the bulk separately using a defined chemical potential, whereas the EqMD method simultaneously models a bulk in equilibrium with the pore. For each GCMC simulation, 40 000 cycles were used, where each cycle consists of a displacement move, an insertion, and a deletion. For particle displacement moves, the probability of success was set to 0.25.
At higher densities, the GCMC method can become less efficient, as insertion and deletion of particles is more difficult. A particular advantage of the BD-NEMD method proposed is that the same system set up can be used to measure adsorbed densities by turning the external forces to zero, removing additional burden of setting up new simulations.
To calculate chemical potentials, for this particular force field one may use directly a molecular based equation of state (EoS), SAFT-γ Mie. 82 The inputs of this EoS are the same as the Mie potential: where A SAFT is the molar Helmholtz free energy of the fluid.
The correspondence between the results of the equation of state and those from molecular simulations has been discussed elsewhere. 83−90 An additional advantage of employing the SAFT-γ Mie EoS alongside the BD-NEMD method is that it allows accessing Λ ij (eq 2), given the condition that the boundaries of the pore are in equilibrium with the bulk regions: where ∂μ/∂z is either calculated using SAFT-γ Mie EoS or from GCMC. Adsorption of Multicomponent Mixtures. Although GCMC and EqMD could be used to calculate adsorbed densities of binary mixtures, given the additional degree of freedom stemming from considering the compositions, the pure component adsorption isotherms of each component are used as input to predict the adsorption of multicomponent mixtures through the ideal adsorption solution theory (IAST). 91,92 In this paper, the IAST method is implemented using the pyIAST package, where the pure component adsorption isotherms (pressure vs adsorbed concentration) are used as inputs. 93 Details are provided in the Supporting Information.
Systems Studied. Four case studies are chosen in this study: (i) Case I: A pure LJ fluid within a slit pore; as a validation of the proposed BD-NEMD method by comparison to transport diffusivities calculated using EqMD. (ii) Case II: Pure methane within a slit pore; as an example of application to realistic fluids. (iii) Case III: A binary mixture of two methane-like fluids, with one having realistic parameters (as in Case II) and the other having augmented interactions with the pore; as an example application on the effect of the solid−fluid interaction on selectivity. (iv) Case IV: A binary LJ mixture, where one species is the same as Case I, and the other has a radius 30% larger; as an example application on the effects of size differences on self-and mutual diffusivities in binary mixtures All parameters for the four systems studied are presented in Table 1. All walls are modeled as FCC lattices, and for all cases, the lattice constant is a 2 2 Wall σ = .
Simulation Details. In all simulations, a system is set up with a solid pore of L z = 12 nm positioned in the middle of a 36 nm simulation box, thus having two bulk regions on each side, each being 12 nm in length. The pore region is an FCC smooth slit pore. Details of its structure and the definition of the pore height are given in the Supporting Information. Lengths of the box (and the pore) in the x and y dimensions are the same, being greater than 10 σ. The large values of pore The Journal of Physical Chemistry B pubs.acs.org/JPCB Article dimensions and simulation box sizes are chosen to minimize finite size effects. In particular, the ratio of particle size to pore length (σ/L) is always less than 0.05, and the ratio of pore height to pore length (H/L) is always less than 1. From previous studies, 94 it is known that these dimensions are far from those resulting in significant finite size effects. The total void volume of the simulation box is known and fixed, and particles are added to a target global fluid density. Equilibrium molecular dynamics simulation are then run in the NVT ensemble. Thus, the system equilibrates, and the adsorption isotherms can be assessed.
With the system equilibrated, the final configuration is used for two purposes. First, the same system is used as the initial configuration for the BD-NEMD simulations. Second, the pore section of the final configuration is isolated (the bulk regions removed and periodic boundary conditions imposed in the xand z-directions) and used in EqMD simulations to calculate Λ ij in an essentially infinite pore setup. The EqMD and NEMD simulations are then used to calculate transport diffusivities which were compared.
Each simulation is run for 10 million time steps of 2 fs. The first 4 ns are used for equilibration, or reaching a steady state, and the remaining 16 ns were analyzed as the main production run.
For both the BD-NEMD and EqMD methods, following previous studies, 37,57 simulations were run in the NVT W ensemble, where T W refers to the temperature coupling of the solid particles only. This means that no temperature coupling is used for the fluid as to avoid influencing its dynamics, and instead energy was added or removed using wall particles as a thermostat, i.e., the excess energy input of the external force is removed by the walls. To implement this, the pore solid particles are allowed to vibrate about their equilibrium position using position restraints of the harmonic form, with a bonding potential of 10 000 kJ mol −1 nm −2 and an equilibrium bond distance of σ wall . By only applying temperature coupling to the solid particles, the temperature of both solid and fluid are kept constant about the equilibrium temperature, which can be found in the Supporting Information. Moreover, to ensure the positions of the fluid were not modified, no center of mass motion removal was used.
The BD-NEMD simulations are run using a modified version of GROMACS/5.1.2 95 and for each simulation at each state point, 20 different simulations were run using applied accelerations in the range of 0−0.004 nm ps −2 to keep the perturbations of the system within the "linear regime". EqMD simulations are run using GROMACS/2018 96 to ensure results are not subject to software bias. All visualizations of simulations have been rendered using the VMD package. 97 ■ RESULTS AND DISCUSSION Case I: Pure LJ Fluid. For the pure LJ system in a slit pore of height 2 nm (6.67 σ), simulations are run from the dilute limit to highly concentrated systems, i.e., ρσ 3 = 0.02−0.85 at T = 1.5 ε/k B , which is above the critical temperature of the LJ fluid (T c = 1.31 ε/k B ). 98 The coexisting fluid density at the supercritical fluid−solid transition for a bulk LJ fluid is at ρσ 3 = 1.015 99 which set the upper limit of densities to 0.85, ensuring the fluid inside the pore did not undergo a phase transition to a solid phase. Figure 2 (right) shows the measured adsorbed concentrations for all bulk concentrations studied. It could be clearly seen that the pore is very selective at lower densities, while at higher densities the concentration is only around five percent less than bulk concentrations. The latter underestimation is presumably an artifact of the particular definition of the pore height used in this work.
Adsorption Isotherms. The BD-NEMD method can be used to measure both transport diffusivities (eq 15) and collective diffusivities, Λ ij , (eq 17). While the first calculation is direct, for the latter, one needs to find the chemical potential gradients across the pore. Otherwise, Λ ij can also be directly calculated using EqMD (eq 2). As a test of robustness of the BD-NEMD method developed, Λ ij is compared with both methods.
Furthermore, chemical potentials are required to evaluate the Darken factors, so Λ ij could be used to calculate D t (eq 13). Chemical potentials are assessed using both GCMC and the SAFT-γ Mie EoS. While using the EoS, the chemical potential is calculated using the bulk density at equilibrium with the adsorbed density, whereas in GCMC (μVT ensemble) chemical potential is an input from which adsorbed densities can be determined. Although the two methods are fundamentally different to each other, a quantitative agreement between them can be seen in Figure 3. The disagreement at higher densities is attributed to sampling inefficiencies in GCMC. The equation of state is computationally much more efficient; thus, the SAFT-γ Mie EoS is chosen as the preferred method of measuring chemical potentials and Darken factors.
The relationship between chemical potential and adsorbed densities can be used to calculate chemical potential gradients across the pore from the BD-NEMD simulations, which can then be used to determine Λ ij using eq 17. Additionally, the Darken correction factor is also calculated from the same relationship: Γ ii is calculated as a function of the adsorbed concentration using SAFT-γ Mie EoS, and the results can be seen in Figure 4 (right). The Darken factor approaches the value of 1 at infinite dilution, corresponding to the expected value for an ideal gas. As more particles are added at lower densities, the presence of attractive forces lead to a decrease in the value of the Darken factor. The Darken factor increases as further insertion of fluid particles becomes less favorable, culminating in very large increases in the chemical potential. The density profile inside the pore can be seen in the Supporting Information. At low The Journal of Physical Chemistry B pubs.acs.org/JPCB Article densities, there is only a small layer of adsorption nearest to the surface. At high density, the pore is saturated, and the fluid is highly ordered. The onset of saturation occurs at ρσ 3 = 0.4, after which particle insertions become less favorable and the chemical potential increases.
EqMD. For each value of the density of fluid inside the slit pore, 60 independent equilibrium simulations, or 300 million time steps in total, were run. Self-diffusivities were calculated using eq 1. Then, the mean square displacement (MSD) of the center of mass of the fluid in the two dimensions (plane) parallel to the pore surface is measured as a function of time. As can be seen in eq 2, in the limit of an infinite time, collective diffusivity, Λ ii is given by the slope of the linear line describing MSD as a function of time. Figure 5 shows the results for one density (ρσ 3 = 0.5), and although the average value of all 60 simulations is linear, each simulation has a drastically different MSD profile relative to the average. Thus, the errors associated with measuring collective diffusivities can be particularly large. This is found to be true across all concentrations.
With the collective diffusivities measured from EqMD and the Darken factors from adsorption isotherms, it is possible to calculate transport diffusivities, D t = ΓΛ (eq 11).
NEMD. For different densities, 20 different forces are applied at the boundary and fluxes and concentration gradients are evaluated. As can be seen in Figure 6a,b, there is a linear relationship between the applied force and the flux as well as the density gradient, i.e., the system is kept within the linear response regime. Frentrup et al. 37 observed a nonlinear response when employing forces an order of magnitude higher than that applied in this work.
In Figure 6a,b, the fluxes and density gradients are shown for different applied forces, respectively. The symbols are colored based on the equilibrium adsorbed densities, i.e., where the force is null. Figure 6c,d shows the contour maps of fluxes and density gradients as functions of the force and the fluid density. Although one would expect that increasing the external force leads to greater flux and density gradients, for a given applied force the value of flux (or density gradient) peaks at intermediate loadings. This could be explained by noting that at higher densities applying a force to particles at the boundary pushes them into a dense fluid region and hinders the flow. However, this does not mean that the transport diffusivity is highest at intermediate concentrations. The value of flux should not be taken as an indication of faster transport. It is only by relating fluxes to density gradients where a true measure of transport is given, which can be seen in Figure 7. The aforementioned figure clearly highlights that increasing pore loading leads to higher transport diffusivities.
Comparison between EqMD and BD-NEMD. A summary of all transport coefficients measured using the EqMD and BD-NEMD methods is presented in Figure 8. From the EqMD method, self-diffusivity (D self ) and collective diffusivity (Λ EqMD ) are calculated using eqs 1 and 2, respectively, and transport diffusivity (D EqMD t ) is calculated using eq 13, i.e., by multiplying Λ EqMD by the Darken factor, Γ. From the BD-NEMD method, Λ BD-NEMD and D BD-NEMD t were calculated using eqs 17 and 15, respectively.
As can be seen, there is quantitative agreement between the BD-NEMD method and the benchmark equilibrium simulations across all densities. The BD-NEMD method captures the same trend as the EqMD simulations; however, both transport and collective diffusivities are underestimated in the nonequilibrium simulation. This underestimation is about 20% in the value of Λ and 15% in the value of D t . The error is greatest at lower densities, which is presumably caused by uncertainties in the equilibrium simulations, as the error bars are very large and the values obtained using the NEMD method all lie within the error associated with EqMD.
It is important to note that while the coefficients approach the same value at infinite dilution, D ρ→0 s = Λ ρ→0 = D ρ→0 t ≈ 3.5 σ(ε/m) 2 , upon increasing densities, these coefficients show different trends. Self-diffusivity decreases with increasing densities and is roughly 20 times less at the highest density than at the infinite dilution limit, as particles in denser phases have smaller velocities and less free paths to diffuse uninterruptedly. This trend is not observed in the collective diffusivities, as the value of Λ peaks at intermediate densities of ρ Ads σ 3 = 0.45 and then slightly decreases at higher density (see Figure 6). Given that the Darken factor substantially increases at higher densities, transport diffusivities show the opposite trend to the self-diffusivities being up to 15 times larger at the highest density relative to the lowest. The fact that transport diffusivities can be orders of magnitude larger than selfdiffusivities stems from the very significant collective motion of the fluids as a consequence of smoothness of the surface of the slit pore, and any individual movement of adsorbed molecules correlates with the movement of all other molecules in the system.
The key difference between this work and previous BD-NEMD approaches in literature is in the way density gradients are defined in the calculations of the transport diffusivities. As mentioned previously, a common practice is to use the unambiguous bulk reservoir densities. However, as can be seen in Figure 9, using bulk densities as the driving force for calculating transport coefficients inside the porous regions leads to drastically different profile for transport diffusivity as a  The Journal of Physical Chemistry B pubs.acs.org/JPCB Article function of density. Bulk density gradients do not take into account the effect of pore adsorption on transport and/or any pore entrance effects that may be present. Using bulk density gradients leads to significant overestimation of transport diffusivity at low loadings and underestimations at intermediate loadings, with a clear minimum at around ρ Ads σ 3 = 0.4. This minimum is neither observed for the BD-NEMD method developed in this work nor the EqMD benchmark (see Figure  9). Case II: Methane. To exemplify how this methodology could be employed to investigate transport of real fluids, the BD-NEMD method is further tested on a system consisting of supercritical methane in a slit pore at T = 300 K. Essentially, the system resembles that of Case I but the fluid force field can be traced back to a realistic model. The Mie parameters describing the intermolecular interactions have been previously optimized to reproduce vapor−liquid equilibrium 100 (see the Supporting Information). The quantitative agreement between results obtained from MD with those obtained from EoS is a unique trait of the SAFT implementation used describing the macroscopic properties of the Mie intermolecular potential, allowing for fast and accurate description of the free energy, and thus chemical potentials of the system.
The top-down coarse-graining technique used to parametrize the force field is effective at producing robust models with transferability and representability, and can be used with confidence to describe transport and adsorption proper-    The slit pore is composed of particles explicitly modeled in an FCC lattice, with self-interaction parameters described using an ad hoc LJ potential, resembling an organic substrate. The height of the pore is 2.6 nm.
The relationship between adsorbed and bulk densities and Darken factors is seen in Figure 10 (left) and in the Supporting Information, where it is seen that at low bulk densities the adsorbed densities are higher than the bulk densities. The pore is slightly wider making surface adsorption less pronounced than the adsorption for Case I. Furthermore, relative to the solid particles used in Case I, the solid particles here are larger with a smaller ε, thus leading to weaker attractive energy density, and thus adsorption capacity.
The comparison between the BD-NEMD and EqMD techniques is presented in Figure 10 (right), showing the remarkable agreement between the two methods, with agreements from the dilute limit to very dense fluids spanning orders of magnitude.
Case III: Transport of a Binary Mixture of Two Methane-Like Fluids with Different Adsorption. Case II was extended by introducing an additional fluid component with similar characteristics to the methane model, henceforth described as MA. This new fluid, MB, has the same intermolecular interactions except that its interaction with the solid pore is less favorable. The cross species energy, ε ij , was 20% less than that of the MA−solid interaction (Table 1). In the bulk, there is no difference between the two fluids, and the mixture can be seen as having the same bulk properties as pure MA at the same total density. However, the reduction in the value of the cross-species interaction results in slightly different adsorption isotherms, and different pore selectivities.
In order to observe the adsorption behavior of the binary mixture, pure component adsorption isotherms are evaluated and can be seen in the Supporting Information. Ostensibly, the adsorption isotherms of the pure fluids look very similar, but given the different interactions with the pore, in a binary mixture it is expected that the pore would be more selective toward species MA. IAST is used to correlate the adsorption isotherms of the mixture using the pure component adsorption isotherm, estimating pore loadings using bulk pressures and compositions. 91,92 The predictions of the IAST were then validated against EqMD simulations of binary mixtures, results that can be found in the Supporting Information. Generally, higher partial pressure of each species leads to higher adsorption of said species.
Pore selectivity toward species MA is defined as where x i is the mole fraction of species i and the superscripts "Ads" and "Bulk" refer to compositions inside the pore and in bulk. A value greater than one refers to a pore that is more selective toward species MA. Figure 11 (left) showcases pore selectivity, with selectivities of up to 16% at the lowest pressures (or bulk densities). With increasing pressure the pore become less selective. The main difference between the adsorption behavior of the two species is the amount adsorbed nearest to the pore surface, where species MA is adsorbed considerably more strongly than species MB, given the less favorable interaction between species MB and the solid. Free energy calculations found in the Supporting Information further showcase the stronger adsorption of species MA. With increasing densities, the first few monolayers near the pore surface become filled, and adsorption occurs in the middle of the pore, which is weakly affected by the solid surface, thus impartially allowing both species to adsorb in the middle resulting in a lower overall selectivity.
The region outlined by the red line in Figure 11 (left) was chosen as the subspace to measure the transport diffusivity matrix of the mixture. This region follows the following constraint:   To measure transport diffusivities from the EqMD method as benchmark cases, Darken factors were calculated using ∂ ln f/∂ρ as described in eq 13 ( Figure 12). Generally, the values of the self-species Darken factor, Γ ii , increases with increasing concentrations of species i. Moreover, Γ ii is much larger that the cross-species Darken factor, Γ ij , indicating that the fugacity of one species is not as strongly correlated with the changes of composition of the other species as it is with changes of compositions of itself. The transport diffusivity matrix is evaluated using both the EqMD and BD-NEMD methods. For the binary mixture, there are four Λ ij and D ij t values to be evaluated with eq 2 used to measure collective diffusivity matrix, [Λ], from the EqMD method. The matrix, along with Γ ij , is used to calculate the transport diffusivity matrix, [D t ] using eq 13. The values of Λ ij at different compositions can be found in the Supporting Information. In general, for this particular system, given the identical interactions between all fluid particles independent of species, it is found that Λ MA,MA ≈ Λ MA,MB and Λ MB,MB ≈ Λ MB,MA , i.e., the flux of each species is equally influenced by chemical potential gradients of either species. This trend is also true for the transport diffusivities measured, and these can be seen in Figure 13, where filled symbols are transport diffusivities measured using the BD-NEMD method, and the empty symbols are those measured from the EqMD method. The first column describes the selftransport diffusivities, and the second column shows the cross coefficients, i.e., the contribution to the flux of species i due to concentrations gradient of species j. . This does not necessarily mean that species MB travels faster than species MA inside the pore, as a concentration gradient in each species causes significant fluxes in the other species. Nevertheless, the high transport coefficients of MB can be explained by the lower selectivity of species MB, as MA particles adsorb more strongly to the surface and are slightly slowed down. If we were interested in evaluating transport selectivity, then transport diffusivities could be used in continuum models to assess changes in composition downstream of a pore.
Commensurate with previous cases, it is seen that transport diffusivities measured using equilibrium methods have large uncertainties. Moreover, the boundary driven method consistently overestimates the diffusivity, with an average of 18% for all elements at all compositions. However, the trends are qualitatively consistent with the values measured using the BD-NEMD method within the uncertainty of the EqMD measurements. This deviation is not significant given that values of transport diffusivities span orders of magnitude with changing concentrations, as previously discussed for Cases I and II.
Case IV: Binary Mixture of LJ Fluid with Size Difference. To understand the effect of size heterogeneity on the elements of the transport diffusivity matrix for a binary mixture, the system studied in Case I was modified by adding another species which is 30% larger than the original LJ fluid, keeping all other self-interaction parameters the same. The pore height is kept the same and the wall particles are thermostated at a temperature of 1.5 ε, resulting in a supercritical fluid. Henceforth the original LJ species will be referred to LJF-1, and the larger species will be referred to as LJF-2. From Table 1, it can be seen that the cross interaction energy, ε ij , between LJF-2 and LJF-1 and between LJF-2 and the solid (LJW) is slightly less than 200 K. This is a consequence of the combining rules used, where the cross interaction energy scales down if the two species have large size differences.
For this system, transport diffusivities were measured for all compositions where ∑ i ρ i Ads σ i 3 = 0.2−0.7. Compared with Case III, there are slight differences in the transport behavior in this system. To highlight the differences, the concentration constraint used to study Case III, i.e., ∑ i ρ i Ads σ i 3 = 0.5, was also investigated. The results can be seen in Figure 14 (left), where it can be seen that the relationship between transport  , and the self-and cross-transport coefficients of LJF-2 are conversely higher than expected. This is in agreement with previous studies, where strong correlation effects lead to the slowing down of a more mobile species and less strongly adsorbed species by the less mobile species. 31 Without a linear relationship between mole fractions and transport coefficients, pure component transport diffusivities cannot be used to approximate the self-transport coefficients in the binary diffusivity matrix as functions of composition. Interestingly, a linear relationship becomes apparent when transport is measured against volume fractions, v or by multiplying densities by σ 3 . When volume fraction of species i, v i , is zero, the self-transport diffusivity of species D ii t is also zero. The self-transport coefficient linearly increases with increasing volume fractions until the point where v i approaches 1, when its value approaches the pure component transport diffusivity.
As with previous cases, in this binary mixture, self-transport diffusivities are orders of magnitude larger than selfdiffusivities. For the systems presented where reduced density is kept constant, self-diffusivities are independent of composition, having a value of 13 × 10 −5 cm 2 s −1 for LJF-1 and 10 × 10 −5 cm 2 s −1 for LJF-2. Again, this is very different to the transport diffusivities seen in Figure 14, emphasizing the fact that self-diffusivities are not adequate parameters to be used in understanding transport in mesoporous materials.
Moreover, for this system the cross coefficients behave differently from the ones studied in Case III. The cross transport coefficient for each species is not similar to the selfterm. For species LJF-1, the cross term D LJF-1,LJF-2 t is larger than the self-term, D LJF-1,LJF-1 t . For the other species, the opposite trend is observed. As previously described, the mutual diffusivity, D ji t , quantifies the influence of concentration gradients of species i on the flux of species j. To understand how the concentration gradient of species i affects the flow of both species (i and j), the ratio D ii t /D ji t was compared with x i /x j at different compositions. If the values of the diffusivity ratio are the same as the composition ratio, then it can be concluded that the flow is fully mutualized and that there is no transport selectivity. This is because each column of the transport diffusivity matrix describes the resultant fluxes emerging from the same concentration gradient, and if the ratio of the elements of each column of the transport diffusivity matrix ratio is the same as the ratio of molar compositions in the pore, then it implies that the fluid is fully mixed, i.e., "toothpaste" or ideal piston flow. An exhaustive discussion is provided in the Supporting Information. Table 2 shows the comparison of ratios and a quantitative agreement is observed.
Although the transport is fully mutualized, composition still affects the overall transport. Increasing the composition of the lighter species leads to a faster flow, which can be quantified by measuring the total fluxes and the total density gradient across the pore (refer to the Supporting Information). The smaller species acts as a diluting agent for the larger species, enhancing overall transport. In this case, the self-transport coefficient of the smaller species is at least 20% larger than that of the larger one in the pure component limit, and the cross coefficient of the smaller species is 7 times larger. Although these metrics would ostensibly allude to a faster transport of species LJF-1, this is not the correct conclusion, as the fluid is fully mutualized, yet it could be seen that for high density fluids inside slit pores, if the adsorption is ideal and the fluid is well-mixed, then the selftransport diffusivity of each component can be estimated using the linear relationship with respect to the volume fractions. Moreover, if the transport is fully mutualized, then the cross terms can be estimated using the ratio of the molar compositions.
The complete picture of transport diffusivity matrix of this system as a function of the adsorbed concentration can be found in the Supporting Information.

■ CONCLUSIONS
This work was carried out to validate the assumptions currently used in the BD-NEMD methods in measuring transport diffusivity of pure components. By benchmarking values of transport diffusivity obtained using the BD-NEMD method against those from the EqMD method, it was highlighted that current implementations of this nonequilibrium method, where it is assumed that concentration gradients can be calculated using bulk reservoir concentrations on either side of the pore, lead to spurious values of transport diffusivity. By relating bulk concentrations on each side to their adsorbed concentrations at equilibrium and using adsorbed concentration gradients, it was shown that the BD-NEMD method can be used to compute accurately transport diffusivities with much lower uncertainty than the EqMD method.
To relate EqMD collective diffusivities to transport diffusivities, Darken correction factors are required. Commonly, these correction factors are assessed using GCMC methods. In this work, a novel method was introduced, Figure 14. Transport diffusivities, D t , evaluated for Case IV LJ mixture at ∑ i ρ i Ads σ i 3 = 0.5. Here "1" and "2" refer to species LJF-1 and LJF-2 respectively. For this plot, total reduced density was maintained at ∑ i ρ i σ i 3 =0.5. (left) Transport diffusivities measured against mole fraction of species LJF-1, x LJF-1 . (right) D t measured against volume fraction of said species, v LJF-1 . The inset plots show the nonlinear relationship of D ii t vs mole fraction and a linear relationship with volume fraction. The results agree well with the ratio of the molar compositions at different volume fractions at ∑ i ρ i Ads σ i 3 = 0.5. This agreement alludes to the fact that the flow is fully mutualized. The subscripts "1" and "2" refer to species LJF-1 and LJF-2 respectively.
The Journal of Physical Chemistry B pubs.acs.org/JPCB Article whereby the Darken factors are evaluated for systems where intermolecular interactions are described using Mie potentials. Using the same system used for the BD-NEMD method at zero force, the equilibrium bulk compositions for each system are used to calculate chemical potentials from the molecular based equation of state, SAFT-γ Mie. A Python code for the evaluation of the SAFT-γ Mie EoS is available. 105 For dense fluids and liquids, it is shown that transport coefficients are orders of magnitude larger than selfdiffusivities. Although there is no formal expectation that these two quantities match except at the ideal gas limit, it is a frequent working assumption which is proven wrong.
In dense binary systems, the influence of the adsorption selectivity on transport is minimal, at least for modest differences in pore attraction. The influence of the size of the particles is, on the contrary, more pronounced. Larger particles dominate the transport, and it is seen that the transport cross coefficients become relevant. These significant values of cross coefficients result in mutualized flow in slit pores at high densities, i.e., no transport selectivity is observed. We do make the caveat that these heuristic observations are relevant only to smooth pores, and a significant difference is seen upon the consideration of transport in rugous nanopores, as we will describe in a future communication.
Full details of the molecular interactions; pore structure and pore height characterization; thermostating; accuracy of the methane model in describing thermodynamic and transport properties; ideal adsorption solution theory for binary mixtures; density profiles inside the pores; enhanced discussion on flow mutualization and transport inside the pores; binary transport coefficients for cases III and IV; free energy calculations of binary mixtures (PDF)